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The description of neuron dynamics can use two distinct representations. On the 
one hand, the membrane potential is the physical variable describing the state of the 
neuron and its evolution is ruled by fundamental laws of physics. On the other hand, a 
neuron is an excitable medium and its activity is manifested by emission of action po- 
tential or "spikes": individual spikes, bursts, spikes trains etc... The first representation 
constitutes the basis of almost all neuron models, and the Hodgkin-Huxley equations 
are, from this point of view, certainly one of the most achieved mathematical represen- 
tation of the neuron [29]. However, neurons communicate by emission of spikes, and 
it is likely that the information is encoded in the neural code, that is, the sequences 
of spikes exchanged by the neurons and their firing times. Since the spikes emission 
results from the dynamics of membrane potentials, the information contained in spikes 
trains is certainly also contained in membrane potential dynamics. But switching from 
membrane potentials to spikes dynamics allows one to focus on information processing 
aspects [26]. However, this change of description is far from being evident, even when 
using simple neuron models (see [38] for a review) . Modeling a spike by a certain shape 
(Dirac peaks or more complex forms), with a certain refractory period, etc .. which 
information have we captured and what have we lost ? These questions are certainly 
too complex to be answered in a general setting (for a remarkable description of spikes 
dynamics and coding see [42]). 

Instead, it can be useful to focus on simplified models of neural networks, where 
the correspondence between the membrane potential dynamics and spiking sequences 
can be written explicitly. This is one of the goals of the present work. We consider 
a simple model of spiking neuron, derived from the leaky integrate and fire model 
[26], but where the time is discretised. To be the best of our knowledge, this model 
has been first introduced by G. Beslon, O. Mazet and H. Soula [50], [51], and we shall 
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call it "the BMS model". Certainly, the simplifications involved, especially the time 
discretisation, raise delicate problems concerning biological interpretations, compared 
to more elaborated models or to biological neurons [14] (see the discussion section). 
But the main interest of the model is its simplicity and the fact that, as shown in the 
present paper, one can establish an explicit one-to-one correspondence between the 
membrane potential dynamics and the dynamics of spikes. Thus, no information is lost 
when switching from one description to the other, even when the spiking sequences 
have a complex structure. Moreover, this correspondence opens up the possibility of 
using tools from dynamical systems theory, ergodic theory, and statistical physics to 
address questions such as: 

— How to measure the information content of a spiking sequence ? 

— What is the effect of synaptic plasticity (Long Term Depression, Long Term Potenti- 
ation, Spike Time Dependent Plasticity, Hebbian learning) on the spiking sequences 
displayed by the neural network ? 

— What is the relation between a presented input and the resulting spiking sequence, 
before and after learning. 

— What is the effect of stochastic perturbations ? Can we relate the dynamics of 
the discrete time BMS model with noise to previous studies on continuous time 
Integrate and Fire neural networks perturbed by a Brownian noise (e.g. [5], [43]) ? 

This paper is the first one of a scries trying to address some of these questions in 
the context of BMS model. The goal the present article, is to pose the mathematical 
framework used for subsequent developments. In section 2 we present the BMS model 
and provide elementary mathematical results on the system dynamics. We show that 
the presence of a sharp threshold for the model definition of neuron firing induces 



4 



singularities responsible for a weak form of initial conditions sensitivity. This effect is 
different from the usual notion of chaos since it arises punctually, whenever a trajectory 
intersects a zero Lebesgue measure set, called the singularity set. Similar effects are en- 
countered in billiards [16] or in Self-Organized Criticality [2], [3], [4]. Applying methods 
from dynamical systems theory we derive rigorous results describing the asymptotic 
dynamics in section 3. Although we show that the dynamics is generically periodic, the 
presence of a singularity set has strong effects. In particular the number of periodic 
orbits and the transients growth exponentially as the distance between the attractor 
and the singularity set tends to zero. This has a strong impact on the numerics and 
there is a dynamical regime numerically indistinguishable from chaos. Moreover, these 
effects become prominent when perturbing the dynamics or when the infinite size limit 
is considered. In this context we discuss the existence of a Markov partition allowing 
to encode symbolically the dynamics with "spike trains" . In section 4 we indeed show 
that there is a one to one correspondence between the membrane potential dynamics 
and the sequences of spiking patterns ("raster plots"). This opens up the possibility 
to use methods from ergodic theory and statistical mechanics (thermodynamic for- 
malism) to analyse spiking sequences. This aspect will be the central topic of another 
paper. As an example, we briefly analyze the CEise of random synapses and inputs on 
the dynamics and compare our analysis to the results obtained by BMS in [51], [50]. We 
exhibit numerically a sharp transition between a neural death regime where all neurons 
are asymptotically silent, and a phase with long transient having the appearance of a 
chaotic dynamics. This transition occurs for example when the variance of the synaptic 
weights increases. A further increase leads to a periodic dynamics with small period. In 
the discussion section we briefly comment some extensions (effect of Brownian noise, 
use of Gibbs measure to characterize the statistics of spikes) that will be developed in 
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forthcoming papers. 

Warning This paper is essentially mathematically oriented (as the title suggests), 
although some extensive parts are devoted to the interpretation and consequences of 
mathematical results for neural networks. Though the proof of theorems and the tech- 
nical parts can be skipped, the non mathematician reader interested in computational 
neurosciences, may nevertheless have difficulties to find what he gains from this study. 
Let us briefly comment this point. There is still a huge distance between the com- 
plexity of the numerous models of neurons or neural networks, and the mathematical 
analysis of their dynamics, though a couple of remarkable results have been obtained 
within the 50 past years (see e.g. [10] and references therein). This has several conse- 
quences and drawbacks. There is a constant temptation to simplify again and again 
the canonical equations for the neuron dynamics (e.g. Hodgkin-Huxley equations) to 
obtain apparently tractable models. A typical example concerns integrate and fire (IF) 
models. The introduction of sharp threshold and instantaneous reset gives a rather 
simple formulation of neuron activity, and, at the level of an isolated neuron, a couple 
of important quantities such as the next time of firing can be computed exactly. The 
IF structure can be extended to conductance based models [44, 14] closer to biological 
neurons. However, there are quite a few rigorous results dealing with the dynamics of 
IF models at the network level. The present paper provides an example of an IF Neural 
Network analysed in a global and rigorous manner. 

The lack of mathematical results concerning the dynamics of neural networks has 
other consequences. There is an extensive use of numerical simulations, which is fine. 
But the present paper shows the limits of numerics in a model where "neurons" have a 
rather simple structure. What is for more elaborated models ? It also warns the reader 
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against the uncontrolled use of terminologies such as "chaos, edge of chaos, complex- 
ity" . In this paper, mathematics allows us to precisely define and analyse mechanisms 
generating initial conditions sensitivity, which are basically presents in all IF neural 
networks, since they are due to the sharp threshold. We also give a precise meaning to 
the "edge of chaos" and actually give a way to locate it. We evidence mechanisms, such 
as the first firing of a neuron after an arbitrary large time, which can basically exist in 
real neural networks, and raise huge difficulties when willing to decide, experimentally 
or numerically, what is the nature of dynamics. Again, what happens for more elabo- 
rated models ? This work is a first step in providing a mathematical setting allowing 
to handle these questions for more elaborated IF neural networks models [14]. 

1 General context. 

1.1 Model definition. 

Fix iV > a positive integer called "the dimension of the neural network" (the number 
of neurons). Let W be an iV x A?^ matrix, called "the matrix of synaptic weights", 
with entries Wij. It defines an oriented and signed graph, called "the neural network 
associated to W" , with vertices i = 1 . . . A'^ called the "neurons" . There is oriented edge 
J — » i whenever Wij =^ 0. Wjj is called "the synaptic weight from neuron j to neuron 
i" . The synaptic weight is called "excitatory" if Wjj > and "inhibitory" if Wij < 0. 

Each vertex (neuron) i is characterized by a real variable Vi called the "membrane 
potential of neuron i" . Fix a positive real number > called the "firing threshold" . 
Let Z be the function Z{x) = x{x > 0) where x is the indicatrix function. Namely, 
Z{x) = 1 whenever x > 6 and Z{x) = otherwise. Z{Vi) is called the "firing state of 
neuron i" . When Z{Vi) = 1 one says that neuron i "fires" and when Z{Vi) = neuron 
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i is "quiescent". Finally, fix 7 e [0, 1[, called the "leak rate". The discrete time and 
synchronous dynamics of the BMS model is given by: 

V(t + l) = F(V(t)), (1) 
where V = {l^j}^^ is the vector of membrane potentials and F = {Fj} with: 

N 

Fi(V)='rVi{l-Z[Vi]) + Y,WijZ[Vj] + If'"; i = l...N. (2) 

i=i 

The variable Jf^* is called "the external current^ applied to neuron i". We shall assume 
in this paper that this current does not depend on time (see however the discussion 
section from an extension of the present results to time dependent external currents) . 
The dynamical system (1) is then autonomous. 
In the following we shall use the quantity 

N 

l!{V) = J2WijZ[Vj]. (3) 
called the "synaptic current" received by neuron i. The "total current" is : 

7i(V) = 7f(V) + 7r* (4) 
Define the firing times of neuron i, for the trajectory^ V, by: 

rf )(V) = inf It > rf -^)(V), V,it) > 6^ (5) 

^ Prom a strict point of view, this is rather a potential. Indeed, this term is divided by a 

capacity C that we have set equal to 1 (see section 1.2 for an interpretation of equation (1)). 

We shall not use this distinction in the present paper. 

Note that, since tlic dynamics is deterministic, it is equivalent to fix the forward trajectory 

or the initial condition V = V(0). 
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where Tj^ = — oo. 

1.2 Interpretation of BMS model as a Neural Network. 

The BMS model is based on the evolution equation for the leaky integrate and fire 

neuron [26] : 

dVi Vi lAt) 

where t = RC is the integration time scale, with R, the membrane resistance and C 
the electric capacitance of the membrane. Ii{t) is the synaptic current (spikes emitted 
by other neurons and transmitted to neuron i via the synapses j — > i) and an external 
current. The equation (6) holds whenever the membrane potential is smaller than a 
threshold 6, usually depending on time (to account for characteristics such as refractory 
period of the neuron). When the membrane potential exceeds the threshold value, the 
neuron "fires" (emission of an action potential or "spike"). The spike shape depends 
on the model. In the present case, the membrane potential is reset instantaneously to 
a value Vreset , corresponding to the value of the membrane potential when the neuron 
is at rest. More elaborated models can be proposed accounting for refractory period, 
spikes shapes, etc ... [26] 

A formal time discretization of (6) (say with an Euler scheme) gives: 

Viit + dt) = Vi{t) (l - ^) + (7) 
Setting dt = 1^ and 7 = 1 — i, wc obtain. 

This can be interpreted as choosing the sampling time scale dt smaller than all charac- 
teristic time scales in the model, with similar effects of refractoriness and synchronization. 
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V^{t + = -lV,(t) + 



(8) 



C 



This discretization imposes that r > 1 in (6), thus 7 e [0, 1[. This equation holds 
whenever Vi{t) < 9. As discussed in e.g. [31] it provides a rough but reahstic ap- 
proximation of biological neurons behaviours. Note that in biological neurons, a spike 
duration is not negligible but has a finite duration (of order 1 ms). 

The firing of neuron i is characterized by: 



where, from now on, we shall consider that C = 1 and that Vreset, the reset potential, 
is equal to 0. Introducing the function Z allows us to write the neuron evolution before 
and after firing in a unique equation (2). Moreover, this apparently naive token provides 
useful insights in terms of symbolic dynamics and interpretation of neural coding. 

Note that the firing is not instantaneous. The membrane potential is maintained at 
a value 9 during the time interval [r-^'^' , r-'*^^ + 1[. Note also that simultaneous firing of 
several neurons can occur. Moreover, a localized excitation may induce a chain reaction 
where ni neurons fire at the next time, inducing the firing of 712 neurons, etc .... Thus, 
a localized input may generate a network reaction on an arbitrary large space scale, in 
a relatively short time scale. The evolution of this propagation phenomenon depends 
However, this requires a more complete discussion, done in a separate paper [14]. See also 
section 5.6. 



and: 



Viir^''^ +1) = Vreset +Ii{r^''^) 



(9) 
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on the synaptic weights and on the membrane potential values of the nodes involved in 
the chain reaction. This effect, reminiscent of the "avalanches" observed in the context 
of self-organized criticality [1], may have an interesting incidence in the neural network 
(1)- 

2 Preliminary results. 

2.1 Phase space M. 

Since 7 < 1 one can restrict the phase space of (1) to a compact set^ M = [Vmim Vmax]^ 
such that F{M) C M where: 



Vmin = niin(0, 



1-7 



j\Wij<0 



(10) 



and: 



Vrr, 



max(0, 



1-7 



j\Wii>0 



ext 



(11) 



where we use the convention "^j^^Wij = 0. Therefore, .<o = (resp. 

SjIW >o ^^'i-j ~ '^^^ weights are positive (rcsp. negative) and "^j^yy. .<q Wij < 
(resp. J2j\w,j>0^ij^^y 

This results is easy to show. Indeed, assume that for all neurons, Vrnin — — ^max- 
Then, the membrane potential of neuron i at the next iteration is 

= ^Vi{i - z{Vi)) + J2 WijZ{Vj) + /r*. 

j 

Therefore, 



Note that in the original version of BMS, Vi > 0. 
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j\Wij<Q 3\Wij>0 



If Vmin < then, 



Vmin = lVmin+ min 
i=l...N 



3\WiS<0 



j\Wij<0 



and if Vmin = 0) then necessarily min 

i=l...N 



Vrr. 



j\Wij<0 



> and V/ > = 



Similarly, if Vmax > then, 



-fVmaxil-Z{Vi))+ y lVy+/r* < 'YVmax+ max 

jW«>o 



E 

jmj>o 



and if Vmax = 0, then necessarily max 

i=l...N 



Vrr, 



j\W,j>0 



ext 



< and Vj < 



Note that the similar bounds hold if /?^* depends on time. 



2.2 Phase space M. 

For each neuron one can decompose the interval J = [Vmin, Vmax] into Tq U Ti with 
To = [Vminj ^1 = [^,Vmax]- If 1^ € Xq the neuron is quiescent, otherwise it fires. 
This splitting induces a partition V oi M, that we call the "natural partition". The 
elements of V have the following form. Call A = {0, 1}^. Let rj — {rji, . . . ,7]^} £ A. 
This is a A'^ dimensional vector with binary components 0, 1. We call such a vector a 
spiking state. Then = |^ Mrj where: 
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Mr, = {V €M\Vi Sir,,} (12) 

Equivalently, V e Mr] ^ Z{Vi) = rji, i = 1 . . . N. Therefore, the partition V 
corresponds to classifying the membrane potential vectors according to their spiking 
state. More precisely, call: 

25(7,) = {»€{!... AT} |r?i = l}, (13) 

and Vir]) the complementary set {i € {1 . . . N} | r;^ = 0}. Then, whatever the mem- 
brane potential V G Adrj the neurons whose index i € ©(»?) will fire at the next 
iteration while the neurons whose index i € T>{r]) will stay quiescent. In particular, the 
synaptic current (3) is fixed by the domain Mri since : 

/nv) = 71(77) = (14) 

whenever V € Mr}- In the same way we shall write Ii{rf) = If{rt) + /f^*. 

V has a simple product structure. Its domains arc hypcrcubos (thus they arc convex) 
where the edges are parallels to the directions Cj (basis vectors of R^). More precisely, 
for each r] € {0, 1}^, 

N 

Mr) = ^hi, (15) 
where W denotes the Cartesian product. 

2.3 Elementary properties of F. 

Some elementary, but essential properties of F, are summarized in the following propo- 
sition. We use the notation 
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JV 

for the cardinality of "Dirf). This is the number of neurons that will fire in the next 
iteration whenever the spiking pattern is t;. 

Proposition 1 Denote by Frj the restriction of F to the domain Mr] ■ Then whatever 

veA, 

1. Ff) is affine and differentiable in the interior of its domain Mr]- 

2. Fr] is a a contraction with coefficient 7(1 — %) in direction i. 

3. Denote by DFrj the Jacobian matrix o/Fjj. Then DFr] has C{r]) zero eigenvalues 
and N — C{r]) eigenvalues 7. 

4- Call Fr]^i the i-th component o/Fjj then 



F{Mr]) = Fr] 



JV 



N 



llFrj^Ir,,) (17) 



i=l 



where Frj^Io) is the interval [iVmin+Y.jLi WijVj+ir^lO+Y.jLi ^ij^j+If'^l 
and Fjj^i(Ti) is the point X^jLi ^ijVj ■ More precisely, ifC{r]) = k, the image 
of Mr] is a N — k dimensional hypercube, with faces parallel to the canonical basis 
vectors ej for all i ^ D(j]) and with a volume 7^"*^ [9 — Vmin]^~'^ ■ 

According to item (1) we call the domains Mr], "domains of continuity" of F. 



Proof By definition, VV € Mr], Fi(V) = jVi{l - m) + E^Li W^jVj + If'- F is 
therefore piecewise affine, with a constant /^(j?) = 'Ylfj^T){r]) ^i^^ fixed by the 

domain Mr]- Moreover Fjj is differentiable on the interior of each domain Mr], with: 



-^=l5ij[l-Vi]. (18) 
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The corresponding Jacobian matrix is thus diagonal, constant in the domain Mrj , 
and its eigenvalues are 7(1 — rji]. Each eigenvalue is therefore if 7?j = 1 (neuron i 
fires) and 7 if 77, = (neuron i is quiescent). Thus, since 7 < 1, Fjj is a contraction in 
each direction i. Once Air] has been fixed, the image of each coordinate Vj is only a 
function of V^. Thus, if V € Mrj = Hili^W' ^^en Frf^i{V) = Frj^iiVi) and t Tf maps 
the hypercube Mrj = riiLi-^»;i onto the hypercube JliLi ^''y,'(-^»/s)- segments 
Trji with rii = are mapped to parallel segments [yVmin + X^^Li ^ijVj + A'^^*i7^ + 
^J—i WijTjj + If^*[ while each segment Irn with % = 1 is mapped to a point. Thus, if 
C(r;) = A; the image of Mr) is a, N — k dimensional hypercube, with faces parallel to the 
canonical basis vectors e,, where i ^ Dij]) and with a volume 7''^"*^ [6* — I4nm]^~'^- 



Finally, we note the following property. The dynamical system (1) can be defined 
on and the contraction property extends to this space. If one considers the J-ball 



□ 



BMiS) = {V € R^\d(V,M) < 6} then : 



F[BMiS)]cBMiS). 



(19) 



The distance d is, for example : 



d(X,X') = max 



i=l...N 



(20) 



natural in the present context according to property 1 (eq. (17)). 



2.4 The singularity set S. 



The set 
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S={V £M, \3i, Vi = 9}, (21) 

is called the singularity set for the map F. F is discontinuous on S. This set has a simple 
structure: this is a finite union of iV — 1 dimensional hyperplanes corresponding to 
faces of the hypercubes M-q. Though iS is a "small" set both from the topological (non 
residual set) and metric (zero Lebesgue measure) point of view, it has an important 
effect on the dynamics. 

Indeed, let us consider the trajectory of a point V € and perturbations with 
an amplitude < e about V. Equivalently, consider the evolution of the e ball B(V, e) 

o o 

under F. If B(V,e) D 5 = then by definition jB(V, e) cMr), some rj, where M-q 
is the interior of the domain Mrj. Thus, by prop. 1(2) F[B(V,e)] C S(F(V),7e). 
More generally, if the images of i3(V, e) under F* never intersect S, then, at time t, 
F*[B(V,e)] C B(F*(V),7*e). Since 7 < 1, there is a contraction of the initial ball, 
and the perturbed trajectories about V become asymptotically indistinguishable from 
the trajectory of V. (Actually, if all neurons have fired after a finite time t then all 
perturbed trajectories collapse onto the trajectory of V after t + 1 iterations). 

On the opposite, assume that there is a time, to such that F*''(B(V,e)) n <S 0. 
By definition, this means that there exists a subset of neurons {ii, . . . ,ik} and V' € 
B(y, e), such that Z{Viito)) ^ 2(V^'(to)), i £ {ii, . . . , ij. Then: 



F,(V{to))-F,{V'{to)) = 

7 [Vi{to){l - Z{Vi{to))) - v:{to){l - Z{Vl{to)))] + Eje{n,...,i,} ^ij [Z{Vj{to)) - Z{V^{to))] 

In this case, the difference between Fj(V(to)) — Fi(V'{to)) is not proportional to 

Vi{to) — Viito) , for i G Moreover, this distance is finite while \Vi{to) — 
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V^{to)\ < e can be arbitrary small. Thus, in this case, the crossing of S by the e-ball 
induces a strong separation effect reminiscent of initial condition sensitivity in chaotic 
dynamical system. But the main difference with chaos is that the present effect occurs 
only when the ball crosses the singularity. (Otherwise the ball is contracted). The 
result is a weak form of initial condition sensitivity and unpredictability occurring also 
in billiards [16] or in models of self-organized criticality [2], [3]. Therefore, S is the only 
source of complexity of the BMS model, and its existence is due to the strict threshold 
in the definition of neuron firing. 

Note that if one replaces the sharp threshold by a smooth one (this amounts to 
replacing an Heaviside function by a sigmoid) then the dynamics become expansive 
in the region where the slope of the regularized threshold is larger than 1. Then, the 
model exhibits chaos in the usual sense (see e.g. [6], [8]). Thus, in some sense, the present 
model can be viewed as a limit of a regular neural network with a sigmoidal transfer 
function. However, when dealing with asymptotic dynamic one has to consider two 
limits {t — > +00 and slope — » +00) that may not commute. 

3 Asymptotic dynamics. 

We now focus on the asymptotic dynamics of (1). 
3.1 The w-limit set. 

Definition 1 (Prom [33,27]) A point y € M is called an oi-limit point for a point 
X € A4 if there exists a sequence of times {ffc}^^, such that x{ti~) — > j/ as tfe — » +00. 
The w-limit set of x, ijj{x), is the set of all w-limit points of x. The w-limit set of M, 
denoted by f2, is the set n = UxgM'^^^)' 
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Equivalently, J? is the set of accumulation points of F*(A^). In the present case, 
since M is closed and invariant, we have O = Ht^o F*(A1). 

The notion of ui limit set is less known and used than the notion of attractor. There 
are several distinct definition of attractor. For example, according to [33]: 

Definition 2 A compact set A £ M is called an attractor for F if there exists a 
neighborhood U of A and a time N > such that F^(W) C U and 

oo 

A=f]F\U). (22) 

t=0 

Note that from equation (19) one may choose for U any open set such that: 

U D Bm{S), yS > 0. (23) 

In our case A and H coincide whenever A ts not empty. However, there are cases 
where the attractor is empty while the lu limit set is not (see example of Fig. 3.3.1 in 
[33], page 128). We shall actually encounter the same situation in section 3.4. For this 
reason we shall mainly use the notion of w-limit set instead of the notion of attractor, 
though we shall see that they coincide except for a non generic set of synaptic weights 
and external currents. 

3.2 Local stable manifolds. 

The stable manifold of V is the set: 

>V'*(V) = {V' |d(F*(V'),F*(V)) -*0 t-^+oo). (24) 

The local stable manifold VV;Qg(V) is the largest connected component of W*(V) 
containing V. It obeys: 
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(25) 



In the present model, if V has a local stable manifold W|(V) of diameter e then: 



Thus, a perturbation of amplitude < e is exponentially damped and the asymptotic 
dynamics of any point belonging to the local stable manifold of V is indistinguishable 
from the evolution of V. 

In BMS model some point may not have a local stable manifold, due to the presence 
of the singularity set. Indeed, if a small ball of size e and center V intersects <S it will 
be cut into several pieces strongly separated by the dynamics. If this happens, V does 
not have a local stable manifold of size e. According to (26) a point Y £ M has a local 
stable manifold of diameter e if : 



where Ug{S) = {V | d(V,5) < 5} is the (5-neighborhood of 5. This means that the 
dynamics contracts the e ball faster than it approaches the singularity set. A condition 
like (27) is useful for measure-theoretic estimations of the set of points having no stable 
manifold via the Borel-Cantelli lemma. 

In the present context, a more direct approach consists in computing: 




(26) 




(27) 



to>0 t>to 



d(V+,5) = inf min \Viit) - e\, 
t>Oi=l...N 



(28) 
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which measures the "distance" between the forward trajectory V"*" = {V(t)}j>o of V 
and S. One has the following: 

Proposition 2 7/d(V~'~,<S) > e > then V has a local stable manifold of diameter e. 

Proof This results directly from proposition 1. Indeed, if d(V~'~,5) > e, the image of 
the e-ball B(V, e) under F*, belong to a unique continuity domain of F, Vt > and F 
is contracting on each domain of continuity. □ 

In the same way, one defines the distance^ of the omega limit set J7 to the singularity 
set (one may also consider the distance to the attracting set whenever A is not empty): 

d{n,S) = ^M^d{V+,S). (29) 

The distance vanishes if and only if J7 fl 5 0. Thus, if d{f2, <S) > e > any point 
of O has a local stable manifold. In this situation, any e- perturbation about V e 12 is 
asymptotically damped. Note however that d{fi,S) can be positive but arbitrary small 
(see section 5.1). 

3.3 Symbolic coding and Markov partition. 

The partition V provides a natural way for encoding the dynamics. Indeed, to each 
forward trajectory V"'' one can associate an infinite sequence of spiking patterns 
rji, . . . ,?7j . . . where rj^ = {vi;t = Z{Vi{t))}^_^. This sequence provides exactly the 
times of firing for each neuron. It contains thus the "neural code" of the BMS model. 

® Note that this is not a proper distance, since one may have d{A, B) = and A ^ B. The 
fact that d( fl, 5) = if and only if J7n5 7^ is true only because both sets are closed. I thank 
one referee for this remark. 
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In fact, this sequence is exactly what biologists call the "raster plot" [26]. On the 
other hand, knowing the spiking sequence and the initial condition V = V(0) one can 
determine V(t) since: 

t-1 t L-l 

Vi{t) = 7* JJ (1 - rh,k) ViiO) + 7 n - Vi;k)Ii{Vn-i), (30) 

fe=0 ra=l k=n 

where Ii{r]^_i) = ^Zf—i ^ijVj;n-i+Ii^* and where we used the convention -y*~" 
%;fe) = 1 if fi = t. (Note that the same equation holds if /f^* depends on time). 

The term 7* Ilfc=o ^ Vi\k) Vi(0) contains the initial condition, but it vanishes 
as soon as ?7j.fc = 1, some k, (which means that the neuron has fired at least once 
between time and t — 1). If the neuron does not fire then this term is asymptotically 
damped. Thus, one can expect that after a sufficiently long time (of order | iog\7)| )' 
system "forgets" its initial condition. Then, knowing the evolution of V(t) should be 
equivalent to knowing the neural code. However, this issue requires a deeper inspection 
using symbolic dynamics techniques and we shall see that the situation is a little bit 
more complex than expected. 

For this, one first defines a transition graph Q^y^ jext^ from the natural partition V. 
This graph depends on the synaptic weights (matrix W) and on the external currents 
(vector I''^*) as well. The vertices of G(wjo^t) are the spiking patterns r] e A = 
{0,1}^. Thus, one associates to each spiking pattern 77 a vertex in ^(wjext). Let 
r},rj' be two vertices of Q^y^jexty Then there is an oriented edge rj if' whenever 
F{Mr]) n M^r 7^ 0- The transition rj —^ rj' is then called legal. Equivalently, a legal 
transition satisfies the compatibility conditions: 
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(a) i € V{r,) n ^ Ejev{r,) + /f"* > ^ 

(6) i e n :d(,7') ^ E.ei^fr,) w^^j- + ^r* < G 

(c) i e n V{r,') ^ + Ejevir,) + ^f"' > ^ 

(d) i € vin) n -Din') ^ ^v^ + E^ei^cr,) + ^r* < ^ 

(recall that ©(t?) is given by eq. (13)). The transition graph depends therefore on the 
coupling matrix W and the external current I''''*. It also depends on the parameters 
7,^ but we shall omit this dependence in the notation. Note that the transitions (a), 
(b) do not depend on the membrane potential. We denote by ■S'^ jext) the set of right 
infinite legal sequences t)^ = {rii, . . . . . .} and by X'^yy jext) the set of bi-infinite 

sequences t) = {. . . j?^, . . . , rj_irjo»7i, ■ ■ ■ , »7t • • •}• 

This coding is particularly useful if there is a one to one correspondence (except for 
a negligible set) between a legal sequence and an orbit of (1). This is not necessarily the 
case due to the presence of the singularity set. However one has this correspondence 
whenever one can construct a finite Markov partition by a suitable refinement oiV. In 
the present context where the dynamics is not expanding and just contracting, a parti- 
tion Q is a Markov partition if its elements satisfy F(Q„)nQ„/ V{Qn) C Qn'- In 
other words, the image of Qn is included in Q„/ whenever the transition n — » n' is legal. 

V is in general not a Markov partition (except if 7 = and maybe for a non 
generic set of Wij , 7?^* values) . This is because the image of a domain usually intersects 
several domains. (In this case the image intersects the singularity set). Prom the neural 
networks point of view this means that it is in general not possible to know what will 
be the spiking pattern at time t + 1 knowing the spiking pattern at time t. There are 
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indeed several possibilities depending on the membrane potential values and not only 
on the firing state of the neurons. The question is however: knowing a sufficiently large 
(but finite) sequence of spiking patterns is it possible, under some circumstances, to 
predict which spiking patterns will come next ? The answer is yes. 

Theorem 1 Assume that d{f2,S) > e > 0. Then: 

1. Call F* the t-th iterate ofF. There is a finite T, depending on d(i7, iS), such that 
T +00 when d{n,S) — > and such that there exists a finite Markov partition 
for F'^. 

2. O is a finite union of stable periodic orbits with a finite period. These orbits are 
encoded by a sequence of finite blocs of spiking patterns, each bloc corresponding to 
a Markov partition element. 

Proof Fix T > 0. Consider the partition p^'^^ whose elements have the form: 

Mr,^...ri^ = Mr,^ n F"^ (Mt,^) n F'^ {Mr],) n . . . n F""^ (Mn^) . (32) 

By construction F^ is continuous and thus is a contraction from the interior of each do- 
main A^jj^^... 77 .j, into >i»7^ , with |F'^(Xrjjj...rj^)| < 'r'^\Mrig...'n^\, where lA^??^...??^ | < 
lA^jjpl and where | | denotes the diameter. Thus there is a finite 



log(e) - logdX^J) 



\og{d{n,S))-\og{\Mr,^\) 



log (7) 

where [] is the integer part, such that 'iJ^ri^^...ri^, |F"^(A^?7^...jj^)| < e < d{n,S). 

(T\ NT N^T 

Then ' has finitely many domains (2 ). Denote them by TTn, n = 1 . . . 2 

Then, |F'^(7r„)| < e,Vn. 

Since F^(i7 n 7r„) C fl F^(7rn) the points belonging to O f) n„ are mapped, 

by F"^, into a subset of f2 of diameter < e. Since d{f2,S) > e > each point in O 
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has a local stable manifold of diameter e. Thus all points of ¥^'"{0 D nn) belong to 
the same stable manifold. Hence all these points converge to the same orbit in O and 
TTn contains at most one point in i?. Since there are finitely many domains ttu, O is 
composed by finitely many points and since the dynamics is deterministic, O is a, finite 
union of stable periodic orbits with a finite period. If 7r„ n 12 = then this domain is, 
by definition, non recurrent and it is mapped into a union of domains tthj. containing 
a point of i7. For all nn containing a point of O, f'^ (irn) n 7r„' 7^ F"^(7rn) C 7r„'. 
Therefore, "P^^^ is a Markov partition for the mapping F^. □ 

Remarks. 

— Structural stability. There is a direct consequence of the previous theorem. As- 
sume that we make a small perturbation of some Wij's or 7f^*'s. This will result 
in slight change of the domains of continuity of V and leads to a perturbed natural 
partition V'. This will also change the w-limit set. Call the perturbed w-limit set 
O'. If d{0,S) > e > then if the perturbation is small enough such that, for any 
orbit in f2, the perturbed and unperturbed orbit have the same sequence of spiking 
patterns, then the set fi and fi' have the same number of fixed points and their 
distance remains small (it vanishes when the amplitude of the perturbation tends 
to zero). This corresponds to a structurally stable situation. On the opposite, when 
increasing continuously the amplitude of the perturbation, there is a moment where 
the perturbed and unperturbed orbit have a different sequence of spiking patterns. 
This corresponds to a bifurcation in the system and the two w-limit sets can be 
drastically different. 

— Maximal period. The number 
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Trf = 2^ 'og(„ , (34) 

gives an upper bound for the number of Markov partition elements, hence for the 
cardinality of Q and for the maximal period. It increases exponentially with the 
system size N and with log(7) and log(d(J7, iS)). (Note that this time is useful 
essentially when d{f2,S) is small (and lower than 1)). Hence, even if the dynamics 
is periodic it can nevertheless be quite a bit complex. 

Theorem 1 opens up the possibility of associating to each orbit in a symbolic 
orbit constituted by a finite sequence of spiking patterns, whenever d(i?,<S) > e > 0. 
This result is generalized in the section 4.1 and its consequence are discussed. 

3.4 Ghost orbits. 

Before proceeding to the characterisation of the w-limit set structure in the general 
case, we have to treat a specific situation, where a neuron takes an arbitrary large time 
to fire. This situation may look strange from a practical point of view, but it has deep 
implications. Indeed, assume that we are in a situation where we cannot bound the first 
time of firing of a neuron. This means that we can observe the dynamics on arbitrary 
long times without being able to predict what will happen later on, because when this 
neuron eventually fire, it may drastically change the evolution. This case is exactly 
related to the chaotic or unpredictable regime of BMS model. Prom a mathematical 
point of view it may induce "bad" properties such as an empty attractor. We shall 
however see that this situation is non generic. 

Definition 3 An orbit V is a ghost orbit if 3i such that: 
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{i}vt>o,Vi{t) <e 

and : 



(u) lim sup Vi {t) 



Examples. 



1. One neuron {N = 1), Wn = 0, Vreset = and /f* = 9(1 - ^) < 9. Take 
Vi{Q) = 0. Then, from eq. (30), Vi{t) = E1=i 7*""-^* = ^(1 - J*) < and 
limf^+oo Vi(t) = e. Therefore the orbit of is a ghost orbit. If Vi(0) > 6 the 
neuron fires and Vi(l) = I^^*. Thus this point is mapped into M = [O, Z'^^*]. 
If < Vi{0) < e then, Vi{t) = 7Vi(0) + ^(1 - 7*) and the neuron fires after a 
finite time, but then it is mapped to Vi = 0. Thus all points of = [O, 7^^*j are 
eventually mapped to and the orbit of is a ghost orbit. In this case O = {0} 
while A is empty (see [33] page 128 for a similar example). 

2. Two neurons with W22 > 0;O < W12 > (1 — f)9; W21 > and where for simplicity 
we assume that Vmin = {Wn > 0) and /f ' = 0. In this case, if 2 fires once, 
it will fire forever. Then the dynamics of 1 is Vi{t + 1) = "yVi^t) + W12, as long 
as Vi{t) < e. Therefore, if Vi{0) < 6, then Vi{t + 1) = 7*+Vi(0) + W12 ^^l!^' 
as long as Vi{t) < 9. The condition Vi{t) < 9 is equivalent to Vi(0) < f{t), with 
fit) = ^ + T^(l - This function is strictly decreasing if 14^12 > (1 - 7)6' 

and f{t) —00 as t — > 00. Thus, for a fixed W12 > (1 — j)9 there is a r = 

ig(i "f^-^) )] 

— log^^^ — (where [ ] is the integer part), such that VO < t < r, there exists 
and interval Jt = [f{t),f{t - l)[e [0,0] such that VVl(O) € Jt, the neuron 1 will 
fire for the first time at time t. When W12 — > 9{1 — 7) from above, r diverges and 
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one can find an initial condition such that the first firing time of 1 is arbitrary large 
(transient case). This generates a ghost orbit. 

One may generalize these examples to arbitrary dimensions. However, the previous 
examples look where very specific since we had to adjust the parameters to a precise 
value, and the ghost orbit can be easily removed by a slight variation of these parame- 
ters. This suggests us that this situation is non generic. We shall prove this in section 
3.5. 

To finish this section let us emphasize that, though "strict" ghost orbits, having 
the limit t — > oo in the definition, are non generic, it may happen that Vi{t) remains 
below the threshold during an arbitrary long (but finite) time before firing. Then, the 
characterization of the asymptotic dynamics may be out of numerical or experimental 
control. 

3.5 Two theorems about the structure of f2. 

The condition d{f2,S) > e > excludes situations where some points accumulate 
on the singularity set. In these situations, the usual behavior is the following. An 
e-ball containing a point V accumulating on S will be cut in several pieces when it 
intersects the singularity set. Then, each of these pieces may intersects S later on, etc... 
At each intersection the dynamics generates distinct orbits and strong separations of 
trajectories. It may happen that the proliferation of orbits born from an e-ball goes on 
forever and there are examples of such dynamical system having a positive (topological) 
entropy even if dynamics is contracting [35]. Also, points accumulating on S do not 
have a local stable manifold. 
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In BMS model the situation is however less complex, due to the reset term 'yVi{l — 
r)i). Indeed, consider the image of an e ball B(V, e) about some point V. Assume that 
the ball intersects several domains of continuity. Then, the action of F generates several 
pieces, as in the usual case. But, the image of B(V, e) CiMrj is a N — C{r]) dimensional 
domain, whose projection in each direction i such that % = 1 is a point. Thus, even 
if B(V, e) intersects the 2^ domains of V, its image will be an union of 2^ pieces all 
but one having a dimension < A'^. This effect limits the proliferation of orbits and the 
complexity of the dynamics and the resulting structure of the w-limit set is relatively 
simple, even if d{n, S) = provided one imposes some additional assumptions. More 
precisely, the following holds. 

Theorem 2 Assume that 3e > and 3T < oo such that, W € Vi e {1 . . . N}, 

1. Either 3t<T such that Vi{t) > 9; 

2. Or 3to = to(V, e) such that Vt > to, Vi{t) < 6 - e 

Then, fl is composed by finitely many periodic orbits with a finite period. 

Note that conditions (1) and (2) are not disjoint. The meaning of these conditions 
is the following. We impose that either a neuron have fired after a finite time (uni- 
formly bounded, i.e. independent of V) or, if it does not fire after a certain time it 
stays bounded below the threshold value (it cannot accumulate on 9). Under these 
assumptions the asymptotic dynamics is periodic and one can predict the evolution 
after observing the system on a finite time horizon T, whatever the initial condition. 
Note however that T can be quite a bit large. 



The proof uses the following lemma. 
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Lemma 1 Fix T a subset of {1 . . . N} and let T he the complementary set of T . Call 

(i) Mi G JP, 3t < T, such that Vi{t) > 6 

(ii) Vj € 3*0 = to(V,i) < oo, such that Wt > to, Vj{t) <9 -e 
then u!{rj7^T,e)! the oj-limit set of rj^^T,€! composed by finitely many periodic orbits 
with a finite period. 



Proof of th. 2 

Note that there are finitely many subsets ^ of {1 . . . A''}. Note also that Pj^ x^e 
/VjT+i,€ and that -fVjT.e C -r:F,T,e' whenever e < e. We have therefore: 



^ ^ U U U ^^.^.^ = [jr^,+oofi- 

:F T>Oe>0 ^ 

But, under hypothesis (1) and (2) of th. 2, there exists e > 0, T < oo such that 
M = \J-prj: T,e where the union on is finite. Since F(A^) C (J^F(ry x^^), Q c 
lj^ci;(ry X j). Under lemma 1 i? is therefore a subset of a finite union of sets containing 
finitely many periodic orbits with a finite period. □ 

Proof of lemma 1 Call Uj: (resp. 11 ^) the projection onto the subspace generated by 
the basis vectors e,, i €: (resp. bj, j € ^) and set = IIjrY (Vf = UfW), 
Fj7 = IIj^F (F^ = n^F). Since each neuron j € ^ is such that: 

Vjit)= I] 7"(^W^,fe2[^fc(t-n-l)l+/f*)<e-e, (35) 

n=0 fe 

for t sufficiently large, (larger than the last (finite) firing time tj), these neurons do 
not act on the other neurons and their membrane potential is only a function of the 
synaptic current generated by the neurons € !F. Thus, the asymptotic dynamics is 
generated by the neurons i € T. Namely, VV e w(r^,7',e), V:F(t + 1) = Fjr[V^(t)] 
and Vjp(t + 1) = F^[Vjf(*)]- One can therefore focus the analysis of the a; limit set 
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to its projection u}j^{rj^^T,e) = ^^<^(^:r,T,e) (and infer the dynamics of the neurons 
j € ^ via eq. (35)). 

Construct now the partition V^'^\ with convex elements given by Mr}Q...r},j, = 
Mt]^ nF"^ [M-q^) nF"2 {-M-r]^) n.. . nF~^ (-^rj^), where T is the same as in the 
definition of rj: ^,^- By construction, F^ is continuous on each element of v'^'^^ and 
fixing Mri„...rij. amounts to fix the affinity constant of F^. By definition of T, DF^|.^, 
the derivative of F^ at V, has all its eigenvalues equal to whenever V e '^j^irj^,T,e) 
(prop. 1.3). Therefore F^[Mr]g...r]^ f^i^j^{rj^,T,e)] is a point. Since 

F^(Xna;^(r^,T,6)) = (\jMrj^...rj^ nujArjr^T,e)) C (J fJ- {Mrj^...Tj^ n ujj.{ry: 

the image of ^^(/V.T.c) under F^ is a finite union of points belonging to M. Since, 
u)j:{rj: X e) is invariant, this is a finite union of points, and thus a finite union of 
periodic orbits with a finite period. The dynamics of neurons € is driven by the 
periodic dynamics of firing neurons and, from eq. (35) it is easy to see that their 
trajectory converges to a constant. □ 

Remark. In the theorem, we have considered the case d(/2, 5) = as well. One sees 
that there is no exponential proliferation of orbits after a finite time corresponding to 
the time where all neurons satisfying property (1) have fired at least once. Indeed, then 
the reset term project a convex domain onto a point, and this point cannot generate 
distinct orbits. As discussed above the effect of S is somehow cancelled by the reset 
intrinsic to BMS model. Note however that there are at most 2^"^ points in i7, and 
this number can be quite a bit large. 

The situation is more complex if one cannot uniformly bound the first time of firing 
as already discussed in section 3.4. Assumptions (1), (2) of theorem 2 leave us on a 
safe ground but are they generic ? Let us now to consider the case where they are not 
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satisfied. Namely Ve > 0, VT < oo, 3V e >1, 3i e {1 . . . N} such that Vt < T, Vi{t) < i 
and Vio, 3t > such that Vi{t) > 9 - e. Call: 



BT,e = < 



(i) yt<T,v,{t)<e 

V e X|3i, such that : } (36) 

{n)yto,3t>to,Vi{t)>e-e. 



We are looking for the set of parameters values (W, I ) such that the set: 



s = n n (37) 

T>Oe>0 

is non empty. Note that Bt+i,^ ^ '^T.e- Thus, B = ne>o B°°t<^- ^^^^ looking for 

points V such that \/t > 0, Vi{t) < 9 and limsup(_>co V^(^) = 9. Therefore, B is exactly 
the set of ghost orbits. 

We now prove that B is generically empty. Actually, we prove a more general result 
namely that d{fi,S) is generically non zero. Before this, we have now to provide a 
definition of "generic" . For this, we shall assume from now on that the synaptic weights 
and inputs belong to some compact space H C '^'^ . This basically means that the 
Wjj's (7f^*'s) are bounded (or have a vanishing probability to become infinite if we deal 
with random matrices/inputs). One can endow Tl with a probability measure having 
a density with respect to the Lebesgue measure. This corresponds to choosing the 
synaptic weights and external currents with some probability distribution, as we shall 
do in section 5.1. We say that a subset J\f Cfi is "non generic in a measure theoretic 
sense" if this set has zero measure. This means that there is a zero probability to pick 
up a point in Af by choosing the synaptic weights and external currents randomly. We 
say that it is "non generic in a topological sense" if it is the complementary set of a 
countable intersection of dense sets [33]. This definition corresponds to the following 
situation. If we find a point belonging to Af then a slight perturbation of this point 
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leads out of A/^, for any perturbation that belongs to an open dense set. In other words 
one can maybe find perturbations that leave the point inside N but they are specific 
and require e.g. precise algebraic relations between the synaptic weights and/or input 
currents. These two notion of genericity usually do not coincide [33]. 

Theorem 3 The subset of parameters (W, P^*) € Ti such that d(Q,S) = is non 
generic in a topological and measure theoretic sense. 

Remark Since this result holds for the two distinct notions of genericity we shall 
use the term "generic" both in a topological and in a measure theoretic sense, without 
further precision in the sequel. 

Proof Take V € J7 such that d(V, S) = 0. Then, there exists i € {1 . . . N} such that 

inft>o \ Vi{t) — 6\ = 0. We shall consider separately two cases. 

1. Either 3B < oo and a sequence {tk}i~>o such that Vi(tfc) = 9 and 5^ < B,Vk > 0, 
where 5k = tfe+i -tfe. 

2. Or V is a ghost orbit. This includes the case where 5^ defined above is not bounded, 
corresponding to having limt_>+tx) Vi{t) = 6, but also the case where Vj(t) has no 
limit, and where lim supj_>_|_;^ ^i(i) = ^ as in definition (3). 

Case 1 According to eq. (30), the condition V^i(ifc+i) = 9 writes: 

Vi{tu+^) = J2 7"Uf (ifc+i - n - 1) + ir') = 6, (38) 

n=0 

since t^ is a firing time. Note that we have used the notation /f (t) instead of the 
notation I?{rif), used in eq. (30), for simplicity. 

The synaptic current I^ takes only finitely many values a^i = X^j£x>(?7,) 
where I is an index enumerating the elements of "P (Z < 2^). Thus, the Ui-i's are only 
functions of the Wij 's and they do not depend on the orbits. One can write: 
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ra=0 Z=l 



where: 



a;i;j(ifc+i) = X] '^"^ [^''(^fc+i - n - 1) = ai;j] , (40) 



n=0 



where x is the indicatrix function. One may view the hst {a;j;;(t/j_|_i)}^_^ as the com- 
ponents of a vector 'x.i{tk+i) £ R -In this setting, relation (38) writes: 



^mitk+i) = - (41) 



z=i 



since I^^* does not depend on time. Equation (41) defines an afiine hyperplane Pj ^ in 

R2 . 

Call Qj^fe the set of a;j;;(tfc_|_i)'s. This is a finite, disconnected set, with #Qi^k = S^*" 
and whose elements are separated by a distance > 7'^''. Moreover, the a;j;;(tfe_|_i)'s are 
positive. For each k they obey: 

E-m(*^+i)=E^" = t3t (''^ 

1=1 n=0 

This defines a simplex and Qi^^ belongs to this simplex. Note does not depend 
on the parameters W, I^^*. However, the set of a;i;i(t/j+i)'s values appearing in eq. (41) 
is in general a subset of ^ depending on (W, I^^*) 

Now, eq. (41) has a solution if and only if P^ /^ fl Qj 7^ 0. Assume that we have 
found a point R = (W, l^''*) in the parameters space H such that Pj^fc fl Qi^k 0i for 
some k. Since Qj^fc is composed by finitely many isolated points, since the a^i's depend 
continuously on the Wij's and since the affine constant of the hyperplane Pj ^ depends 
continuously of /f^*, one can render the intersection Pij^ <^Qi,k empty by a generic (in 



33 



both sense) small variation of the parameters Wij,If . Therefore, the sets of points 
in 7i such that Pj ^ fl Qj ^ ^ 0, for some k, is non generic. Since we have assumed 
that the S^'s are uniformly bounded by a constant B < oo, the condition 3k such that 
Vi{ti^) = corresponds to a finite union of non generic sets, and it is therefore non 
generic. 

Note that if 5k is not bounded then the set of values x^i = X^^o '>'"x (^fe+i — n — 1) = a^i^ 
takes uncountably many values. If 7 is sufficiently small this is Cantor set and one can 
still use the same kind of argument as above. On the other hand, if 7 is large this 
set fills continuously the simplex Xi-i = and one cannot directly use the 

argument above. More precisely one must use in addition some specificity of the BMS 
dynamics. This case is however a sub case of ghost orbits. Therefore we treat it in the 
next item. 

Case 2. We now prove that ghost orbits are non generic. For this, we prove that 
if = (W, I) is a point in Ti such that the set B defined by eq. (36) is non empty, 
a small, generic, perturbation of R leads to a point such that B is empty. Thus, B is 
generically empty in both sense. 

Fix e and take V € Boo,e (def. (36)). Then there is a to such that 6-e< Vi{to) < 6. 
Without loss of generality (by changing the time origin) one may take to = 0. Then, 
firom eq. (30), Wt > 0, 

t t 
l\e - e) + ^ 7*""^i(ri - 1) < Vi{t + 1) < 7*e + XI 7*""-?i(n - 1), 

n=l n=\ 

where we have set 7j(n — 1) = Ii{r]n-i) to shorten the notations. Thus, Vi{t) belongs 
to an interval of diameter 7*e. Since e can be arbitrarily small, and t arbitrarily large 
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we have only to consider the orbits such that Vi{0) = 9, for some i. There are finitely 
many such orbits. 

Assume that R = (W, I) is such that B is non empty. Then, for some i, Ve > 0, 
there exists to such that: 

to 

6>-e<7*«6i + ^7*"-"7i(n-l) <e, (43) 

n=l 

and Vt > 0, 

t 

^7*-"/i(n-l)<^(l-7*)- (44) 

Tl=l 

Assume for the moment that there is only one neuron i such that infj>o |Vi(t) — ^| = 
0. That is, all other neurons j =^ i are such that Vj (t) stays at a positive distance from 
9. In this case, a small perturbation of the W^j's, where k — 1 . . . N but j 7^ i, or a small 
perturbation of the /|^*'s will not change the values of the quantities rij{t) = Z{Vj{t)), 
t = . . . + 00. In this case, the current /j(n — 1) in eq. (43,44) does not change Vn > 0. 
Therefore there is a whole set of perturbations that do not remove the ghost orbit^ . But 
they are non generic since a generic perturbation involves a variation of all synaptic 
weights Wkj including j = i and all currents as well. 

Now, a small perturbation of some Whi or 7?^* has the following effects. Call V^{t) 
the perturbed value of the membrane potential at time t. 

1. Either Vi > 0, V^'(t) < 9 — eo, for some eo > 0. In this case, condition (43) is violated 

and this perturbation has removed the ghost orbit. Now, since i is not firing, it does 

not act on the other neurons and wc arc done. 
^ For example, there may exist submanifolds in Ti. eorresponding to systems with ghost 
orbits. A possible illustration of this is given in fig. 1, section 5.1 where the sharp transition 
from a large distance d(n,S) to very small distance d{Q,S) corresponds to a critical line in 
the parameters space 7,(7 (see section 5.1 for details). 
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2. Or there is some to such that V^'(to) > The condition (43) is violated and this 
perturbation also removes the ghost orbit. But, neuron i is now firing and we have 
to consider its effects on the other neurons. Note that the induced effects on neurons 
j ^ i is not small since neuron j feels now, at each times where neuron i fires, an 
additional term Wji which can be large. Thus, in this small perturbation 

induces drastic changes by "avalanches" effects. 
Again, we have to consider two cases. 

(a) Either the new dynamical system resulting from this perturbation has no ghost 
orbits and we are done. 

(b) Or, there is another neuron ii {i\ ^ i) having a ghost orbit obeying conditions 
(43,44). But then one can remove this new ghost orbit by a new perturbation. 
Indeed, as argued above, the fact that i is now firing corresponds to adding a 
term Wji to the synaptic current Ij each time neuron i fires. Then, to still have 
a ghost orbit for j one needs specific algebraic relations between the synaptic 
weights and currents which corresponds to a set of parameters of codimension 
lower than 1. The key point is that, following this argument, one can find a 
family of generic perturbation that destroy the ghost orbits of ii without creat- 
ing again a ghost orbit for i. Then by a finite sequence of generic perturbations 
one can find a point in H such that B is empty. 

Finally, we have to treat the case where more than one neuron are such that 
infoo \Vj{t) — 9\=Q. However these neurons correspond to case 1 or to case 2 and one 
can lead them to a positive distance from <S by a finite sequence of generic perturba- 
tions. □ 
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3.6 General structure of the asymptotic dynamics. 

We are now able to fully characterize the u! limit set of M. 

1. Neural death. Assume that /f^* < {1-^)6 and consider the set A^o = {V \Vi < 6, Vi} 
corresponding to states where all neurons are quiescent. Under this assumption on 

If , Mo is an absorbing domain (F(A^o) C Mo) and F (A^o) as t ^ oo. 

Thus, all neurons in this domain are in a "neural death" state in the sense that 
they never fire. More generally, let Mr) be a domain such that 3t > such that 
F^(Mr]) C Mo then all states in M-q converge asymptotically to neural death (un- 
der the assumption If^^ < (1 — 1)9)- Now, if Ut>o^ -'^ ^^"^^ state 
VV G M converges to neural death. Such a condition is fulfilled if the total current 
is not sufficient to maintain a permanent neural activity. This corresponds to the 
previous condition on 7f^* but also to a condition on the synaptic weights Wij. 
For example, an obvious, sufficient condition to have neural death is Vmax < 0- 
More generally, we shall see in section 5.1, where random synapses are considered, 
that there is a sharp transition from neural death to complex activity when the 
weights have sufficiently large values (determined, in the example of section 5.1 by 
the variance of their probability distribution). 

2. Full activity. On the opposite, consider now the domain Mi = {V \Vi > 6, \fi} 
corresponding to states where all neurons are firing. Then, if Vi, ^J—i + 
jext y ^jjjg domain is mapped into itself by F (where F(A^i) is the point 
Y^—x X^jLi + ^''^) all neuron fire at each time step, forever. More gen- 
erally, if Ut>o -^ ^^-^i) then all state W € M converges to this state of 
maximal activity. Such a condition is for example fulfilled if the total current is too 
strong. 
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These two situations are extremal cases that can be reached by tuning the total 
current. In between, the dynamics is quite a bit richer. One can actually distinguish 3 

typical situations described by the following theorem, which is a corollary of Prop. 1, 
th. 1, 2 and previous examples. 

Theorem 4 Let 



V+ = max , (45) 



where: 



1/+= sup limsupVi(i), (46) 
vex t^oo 

he the maximal membrane potential that the neurons can have in the asymptotics. Then, 

1. Either < 6. Then = maxj d{n,S) = 6 — and fi is reduced to a 
fixed point € A^o- [Neural death]. 

2. Or d{n,S) > e> andV+ > 9. Then Q is a finite union of stable periodic orbits 
with a finite period [Stable periodic regime.]. 

3. Or d{f2, S) = 0. Then necessarily > 9. In this case the system exhibits a weak 
form of initial conditions sensitivity, fi may contain ghost orbits but this case is non 
generic. Generically, the ui-limit set is a finite union of periodic orbit. [Unstable 
periodic regime.]. 

Remark 

It results from these theorems that the BMS model is an automaton; namely, 
the value of 77 at time t can be written as a deterministic function of the past spiking 
sequences r]{t — l),r}{t — 2) etc .... However, the number of spiking patterns determining 
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the actual value of r] can be arbitrary large and even infinite, when d{f2,S) = 0. 
Moreover, the dynamics is nevertheless far from being trivial, even in the simplest case 
7 = (see section 5.1). 

4 Coding dynamics with spiking sequences. 

In this section we switch from the dynamics description in terms of orbit to a description 
in terms of spiking patterns. For this we first establish a relation between the values 
that the membrane potentials have on f2 and an infinite spiking patterns sequence, 
using the notion of global orbit introduced in [17]. 

4.1 Global orbits. 

In (30), we have implicitly fixed the initial time at t = 0. One can also fix it at t = s 
then take the limit s — > — oo. This allows us to remove the transients. This leads to: 

Vi{t) = ^n,in,th"l!it-n-l) (47) 

n=0 

where: 

n 

^i{n,t) = Y[{l-m;t-k-i), (48) 

fe=0 

Definition 4 An orbit is global if there exists a legal sequence f) = {»7t}tgiz £ 
I7(yy jext-j such that Vt > 0, Vj(t) is given by (47). 

Remarks 
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1. In (47) one considers sequences r},.t-i~-i where t — k — 1 can be negative, i.e. 
{^*}teiz ^ -^(W,!"^*)- Thus a global orbit is such that its backward trajectory stays 
in M, yt < 0. 

2. The quantity 7rj(n, t) € {0, 1}, and is equal to 1 if and only if neuron i, at time t, 

(k) 

has not fired since time time t — n — 1. Thus, if r^ ' is the last firing time, then 

(k) 

Vi{t) = ~^l"-l!{t - n - 1),T^''^ < t < t('=+^^ is a a sum with a finite 

number of terms. The form (47) is a series only when the neuron didn't fire in the 
(infinite) past. 

Denote by Q the set of global orbits. The next theorem is an (almost) direct trans- 
position of proposition 5.2 proved by Countinho et al. in [17]. However, the paper [17] 
deals with a different model and slight adaptations of the proof have to be made. The 
main difi'erence is the fact that, contrarily to their model, it is not true that every 
point in has a uniformly bounded number of pre-images. This is because F typ- 
ically project a domain onto a domain of lower dimension in all directions where a 
neuron fires (and this effect is not equivalent to setting a = in [17]). Therefore, to 
apply Countinho et al. proof we have to exclude the case where a point has infinitely 
many pre-images. But it is easy to see that in the generic situation of th. 2 any point 
of Q has a finite number of pre-images in Q (since fi has finitely many points). 

The version of Countinho et al. theorem for the BMS model is therefore. 

Theorem 5 . f2 = G for a generic set o/(>V,I^''*) values. 

Remark For technical reasons we shall consider the attractor A definition (eq. 22) 
instead of the w-limit set. But these two notions coincide whenever there is no ghost 
orbit (generic case). 
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Proof The inclusion ^ C ^ is proved as follows. Let V € Q and V = {V(t)}jg|2 be the 
corresponding global orbit. Since, Vt, n, 

N N 

mm ^Wij < Ii{t-n-l) < max^Wij, 

one has 

oo / N \ oo / N 

n=0 y ' j=l ) n=0 y * j=l 

^TOTO =< Vi{t) < Vmax- 

Therefore, V{t) € M C Vt < 0, 5 > 0. Hence V e (X^o^^i^MiS)) and 

e C nt=oF*(^^W)- f^om (19)' nt=oF*(B^W) C A, and C A 



The reverse inclusion ^ C 5 is a direct consequence of the fact that any point of A 
has a pre-image in A. Therefore, VV € A, one can construct an orbit {V(t)}«Q such 
that V(0) = V, V{t + 1) = F(V(t)) and V{t) € A, Vt < 1. This (backward) orbit 
belong to M and the value of V{t) is given by (47). Thus V e Q, so Ac Q. □ 

Remark. Theorem 5 states that each point in the attractor is generically encoded 
by a legal sequence fj. This is one of the key results of this contribution. Indeed, 
as discussed in the introduction, the "physical" or "natural" quantity for the neural 
network is the membrane potential. However, it is also admitted in the neural network 
community that the information transported by the neurons dynamics is contained in 
the sequence of spikes emitted by each neurons. In the BMS model such a sequence 
is exactly given by fi since on the i-th line %;i one can read the sequence of spikes 
(and the firing times) emitted by i. The theorem establishes that, in the BMS model, 
it is equivalent to consider the membrane potentials or the spiking sequences: the 
correspondence is one to one. This suggests a "change of paradigm" where one switches 
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from the dynamics of membrane potential (eq. 1) to the dynamics of spiking patterns 
sequences. This is the point of view developed in this series of papers, where some 
important consequences are inferred. 

5 Discussion 

5.1 Random synapses. 

In this paper we have established general results on the BMS model dynamics, and we 
have established theorems holding either for all possible values of the Wij's and /f^*'s 
or for a generic set. However, and obviously, the dynamics exhibited by the system (1) 
depend on the matrix W (and the input I^^*) and quantities such as d{f2,S) or V"^ 
in th. 4 are dependent on these parameters. A continuous variation of some Wij or 
some ^* will induce quantitative changes in the dynamics (for example it will reduce 
the period or the number or periodic orbits). It is therefore interesting to figure out 
what are the regions in the parameters spaces W, I^^* where the dynamics exhibits a 
different quantitative behaviour. 

A possible way to explore this aspect is to choose W (and/or l'^'^*) randomly, 
with some probability Vyy ('Pjext) having a density. A natural starting point is the 
use of Gaussian independent, identically distributed variables, where one varies the 
statistical parameters (mean and variance). Doing these variations, one performs sort 
of a fuzzy sampling of the parameters space, and one somehow expects the behaviour 
observed for a given value of the statistical parameters to be characteristic of the 
region of W, I^^* that the probabilities Vw,Visxt weight (more precisely, one expects 
to observe a "prevalent" behaviour in the sense of Hunt & al. [30]). 
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Imposing such a probability distribution has several consequences. First, the synap- 
tic currents and the membrane potentials become random variables whose law is in- 
duced by the distribution Pyy lext = Pw^l'"'* ^-iid this law can be somehow determined 
[15]. But, this has another, more subtle effect. Consider the set Sy{ of all possible se- 
quences on A = {1...7V}. Among them, the dynamics (1) selects a subset of legal 
sequences, X'jyy jext), defined by the compatibility conditions (31) and the transition 
graph C/^yy jext). Thus, changing W (I''^*) has the effect of changing the set of legal 
transitions that the dynamics selects. Prom a practical point of view, this simply means 
that the typical raster plots observed in the asymptotic dynamics depend on the Wij 's 
and on the external current I^'^*. This remark is somewhat evident. However, a ques- 
tion is how the statistical parameters of the distribution Pyy je^t acts on the dynamics 
typically observed in the asymptotics (e.g. how it acts on the parameters V"^, d{Q, S)). 
This question can be addressed by combining the dynamical system approach of the 
present paper, probabilistic methods and mean-field approaches from statistical physics 
(see [10,48] for an example of such combination applied to neural networks). A detailed 
description of this aspect would increase consequently the size of the paper, so this will 
be developed in a separate work [15]. Instead, we would like to briefly comment results 
obtained by BMS. 

Indeed, the influence of the statistical parameters of the probability distribution of 
synapses on the dynamics has been investigated by BMS, using a different approach 
than ours. They have considered the case where the M^y's are Gaussian with zero 
mean and a variance cr^, and where the external current was zero. By using a mean- 
field approach they were able to obtain analytically a (non rigorous) self-consistent 
equation (mean-field equation) for the probability xt that a neuron fires at a given 
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time. This equation always exhibits the locally stable solution a; = corresponding 
to the "neural death". For sufficiently large a another stable solution appears by a 
saddle-node bifurcation, corresponding to a non zero probability of firing. In this case, 
one has two stable coexisting regimes (neural death and non zero probability of firing) , 
and one reaches one regime or the other according to the initial probability of firing. 
Basically, if the initial level of firing is high enough, the network is able to maintain 
a regime with a neuronal activity. This situation appears for a sufficiently large value 
of cT, corresponding to a critical line in the plane 7, cr. The analytical form of this 
critical line was not given by BMS. Moreover, the mean-field approach gives information 
about the average behavior of an ensemble of neuraJ networks in the limit N 00. 
The convergence involved in this limit is weak convergence (instead of almost-sure 
convergence). Therefore, it does not tell us what will be the typical behaviour of one 
infinite sized neural network. Finally, the mean-field approach does not allow to describe 
the typical dynamics of a finite sized network. 

To study the finite size dynamics BMS used numerics and gave evidence of three 
regimes. 



— Neural death. After a finite time the neurons stop to fire. 

— Periodic regime. This regime occurs when a is large enough. 

— "Chaos". Moreover, BMS exhibit an intermediate regime, between neural death 
and periodic regime, that they associate to a chaotic activity. In particular, nu- 
merical computations with the Eckmann-Ruelle algorithm [22] exhibit a positive 
Lyapunov exponent. This exponent decreases to zero when a increases, and be- 
comes negative in the periodic regime. 
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Their conclusion concerning the existence of a chaotic regime is in contradiction 
with theorem 4. We would like now to briefly comment this contradiction (a more 
detailed investigation will be done in [15]). The fig. la,b presents the results of a nu- 
merical simulation computing the average distance d{0, S) as a function of 7 and of the 
variance of the synaptic weights. More precisely, we have considered, as BMS, the case 
of Gaussian independent, identically distributed random Jjj's, with zero expectation 
and variance cr^ = (We have adopted the standard scaling of the variance with 
Indeed, in the present case the neural network is almost surely fully connected and the 
scaling ^ is used in order that the probability of the total currents /j has a variance 
independent of N). 

Clearly, the average distance becomes very small when C crosses a critical line in the 
plane C, 7. However, in the numerical experiments of Fig. 1 the smaller measured value 
for the distance is ~ 10~^ for Fig. lb, corresponding to a very large characteristic time 
well beyond the transients usually considered in the numerics (eq. (34). Moreover, the 
average distance approaches zero rapidly as N growths. Thus, there is sharp transition 
from neural death to chaotic activity in the limit iV — > 00, when crossing a critical line 
in the plane C, 7 ( "edge of chaos" ) . This line can be determined by mean-field methods 
analogous to those used in [8] and corresponds to the transition found by BMS [15]. 
In fig. la,b, one also remarks that after the transition d{n,S) growths slowly when C 
increases. For the illustration of this aspect we have drawn the log of the distance in 
fig. la,b. 

Hence, for finite size N the situation is the following. Start from a small variance 
parameter C and increase it, and consider the stationary regime typically observed. 
There is first a neural death regime. After this, there is a regime where the dynamics 
has a large number of periodic orbits and very long transients. This regime is numeri- 
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cally indistinguishable from chaos . In particular, usual numerical methods, computing 
Lyapunov exponents by studying the behaviour of a small ball of perturbed trajectories 
centered around a mother trajectory, will find a positive exponent. Indeed, if the size 
ri of this ball is larger than the distance d{f2, <S) one will observe an effective expansion 
and initial condition sensitivity, as argued in the section 2.4. This will result in the 
measurement of an effective positive Lyapunov exponent, stable with respect to small 
variation of r], as long as >> d{f2,S). Though this exponent is, strictly speaking, 
spurious, it captures the most salient feature of the model: sensitivity to perturbations 
with a finite amplitude. When C increases further, the distance to the singularity set 
increases. There is then a C such that the typical periodic orbit length becomes of the 
order of magnitude of the time range used in the numerical simulation, and one is able 
to see that dynamics is periodic. 

In the light of this analysis we claim that BMS results are essentially correct though 
we have shown that there is no strictly speaking chaotic regime. Moreover, they are, in 
some sense, more relevant than theorems 3,4 as far as numerics and practical aspects 
are concerned. However, the analysis of the present paper permits to have a detailed 
description of the typical dynamics of a given finite sized network (without averaging), 
based on rigorous results. This is useful when dealing with synaptic plasticity and 

^ Moreover, it is likely that the phase space structure has some analogies with spin- 

glasscs[15]. For example, if 7 = the dynamics is essentially equivalent to the Kauffman's 
cellular automaton [32]. It has been shown by Dcrrida and coworkers [20], [21] that the Kauff- 
man's model has a structure similar to the Slierrington-Kirckpatrick spin-glass modcl[39,46]. 
The situation is even more complex when 7 ^ 0. It is likely that wc have in fact a situation 
very similar to discrete time neural networks with firing rates where a similar analogy has been 
exhibited [7], [8]. 
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learning effects where a given pattern is learned in a given network. (This aspect is 
shortly discussed below and will be developed elsewhere). 

5.2 Adding noise to the dynamics. 

It is usual in neural network modeling to add Brownian noise to the deterministic 
dynamics. This noise accounts for different effects such as the diffusion of neurotrans- 
mitters involved in the synaptic transmission, the degrees of freedom neglected by the 
model, external perturbations, etc ... Though it is not evident that the "real noise" is 
Brownian, using this kind of perturbations has the advantage of providing a tractable 
model where standard theorems in the theory of stochastic processes [24] or methods 
in non equilibrium statistical physics (e.g. Fokker-Planck equations [5]) can be applied. 

The addition of this type of noise to the dynamics of BMS model will result, in 
the region where d{n,S) is small, in an effective initial condition sensitivity and an 
effective positive Lyapunov exponent. 

More precisely, consider a noisy version of (1). 

TV 

Vi{t + 1) = ^Vi{t) (1 - Z[Vim + Yl + ^t"\t) + Bi{t); i = 1 ... AT. 

(49) 

where B {B^{t)'\^J.'^^_Q is a Gaussian random process with zero mean and a co- 
variance Cov{Bi{t),Bj{s)) = a^St^s^ij ■ The probability distribution of the stochastic 
process V, on a finite time horizon T, for a fixed realisation of the W can be obtained 
by using a discrete time version of Girsanov theorem [28], [48]. Prom this, it is possible 
to estimate the probability that a trajectory approaches the singularity set S within a 
finite time T and a distance d by using Freidlin-Wentsel estimates [25]. Also, eq. (27) is 
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useful to estimate the measure of points having a local stable manifold. In this context 
one can compute the probability to approach the singularity set within a distance e; 
also one can construct a Markov chain for the transition between the attraction basin 
of the periodic orbits of the unperturbed dynamics. This will be done in a forthcoming 
paper. 

5.3 Time dependent input. 

One may also wonder what happens to the present analysis when a deterministic, time 
dependent external input, is imposed upon the dynamics (the case of a stochastic input 
is covered by eq. (49) above). Away from the singularity set (rf(i?,iS) large) the effect 
of a time dependent input with a small amplitude (lower than d{0,S)) will not be 
different from the case studied in the present paper. This is basically because a small 
input may be viewed as a perturbation of the trajectory, and the contraction properties 
of the dynamics will damp the perturbation as long as the trajectory stays away from 
the singularity set. 

The situation is different if, at some place, the action of the time dependent input 
leads to a crossing of the singularity. This crossing can basically occur with a time 
independent input, but in the time dependent case there is a particularly salient effect, 
that may be easily revealed with periodic external currents. That is resonance effects. 
If the unperturbed trajectory has some typical recurrent time to come close to the 
singularity set, and if the time dependent perturbation is not synchronized with this 
recurrence time, one expects that the contraction effect will damp the perturbation with 
no clear cut "emergent" effect. On the other hand, if the period of the periodic signal 
is a multiple of the recurrence time, there may be a major effect. The result would be a 
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frequency dependent response of the system exhibiting sharp peaks (resonance). This 
statement is actually more than a conjecture. Such resonances effects have indeed been 
exhibited in a recurrent discrete time neural network with firing rates [11], [12], [13]. 
ft has been shown that applying a periodic input is a way to handle the interwoven 
effects of non linear dynamics and synaptic topology. Similar effects should be observed 
in BMS model. 

5.4 Learning and synaptic plasticity. 

What would be the effect of a synaptic weight variations (synaptic plasticity, LTD, 
LTP, STDP, Hebbian learning) on the dynamical system (1) ? These variations corre- 
sponds to moving the point corresponding to the dynamical system in the parameters 
space (W, I**''*). This motion is neither random nor arbitrary. Indeed, assume that 
one imposes to the neural network an input /stimulus I'^''* = {/j?^*(f)}. If^* modifies 
directly the level of activity of neuron i, and acts indirectly on other neurons (provided 
that the synaptic graph is connected) . A simple stimulus can therefore strongly modify 
the dynamics, the attracting set, the distance d{f2,S), etc .... 

In the case where I^^* does not depend on time, the following result follows directly 
from the analysis presented in this paper. 

Theorem 6 For a generic set of values of (W,!^^^), there exists a finite partition of 
M = [jT>n, such that VV € Vn the uj-limit set ofV, a;(V) is a stable periodic orbit, 
with a finite period. This orbit depends on I^^^ . 

Proof Q is generically a finite union of periodic orbits with a finite period. Each of 
orbit n has an attraction basin Vn and the attraction basins consitute a partition of 
M. □ 



49 



This orbit (resp. its coding) may be viewed as the dynamical response of the neural 
network to input I^^*, whenever the initial conditions are chosen € "Dn- In this way, 
the neural network associates to an input a dynamical pattern encoded in the spiking 
sequence of this periodic orbit. In the same way one can associate to a series of inputs a 
series of periodic orbits (resp. codes), each orbit being specifically related to an input. 
This property results directly from th. (6) without particular assumption on the Wij 's. 

However, there might exist a large number of domains Vn and a large number of 
possible responses (orbits). Moreover, an orbit can be complex, with a very long period. 
This is particularly true at the "edge of chaos" . Indeed, consider the case where the 
distance d{n, S) is small, when the input is present. Then, dynamics is indistinguishable 
from chaos and the dynamical "signature" of the input is a very complex orbit, requiring 
a very long time to be identified. In other words, if one imagine a layered structure 
where the present neural network acts as a retina and where another neural network 
is intended to identify the orbit and "recognize" the input, the integration time of the 
retina will be very long at the edge of chaos. On the opposite, one may expect that 
a learning phase allows this system to associate the input to an orbit with a simple 
structure (small period) allowing a fast identification of the input. 

It has been shown, in the case of recurrent neural networks with a sigmoidal transfer 
function [18], that Hebbian learning leads to a reduction of chaos towards a less com- 
plex dynamics, permitting to associate a pattern to simple orbits. The same effect has 
been observed by BMS [51] applying an STDP like rule to the model (1). In both cases, 
it has been observed that a synaptic evolution (Hebb or STDP) leads to associate to 
the input a sequence of orbits whose complexity decreases during the synaptic weights 
evolution. In the present context, this suggests that d{(jj{y),S) increases during this 
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evolution (note that the evolution is entirely dependent on the initial condition, V. ). 

A related question is: how do the statistical properties of raster plots evolve during 
synaptic weights evolution ? This question, and more generally the effect of synaptic 
evolution on dynamics can be addressed using tools from dynamical systems theory, in 
the spirit of the present paper. This will be the subject of a forthcoming paper. However, 
in the next section we mention briefly how tools from ergodic theory (thermodynamic 
formalism) can be used. 

5.5 Statistical properties of orbits. 

As we saw, the dynamics of (1) is a rather complex and can be, from an experimental 
point of view, indistinguishable from chaos. Consequently, the study of the finite evo- 
lution of the membrane potential (resp. the spiking patterns sequence) does not tell us 
what will be the further evolution, whenever the time of observation is smaller than 
the characteristic time of eq. (34). In this sense, the system is producing entropy 
on a finite time horizon. Thus, provided that d{Q, S) is sufficiently small, one can do 
"as if" the system were chaotic and use the tools for analysis of chaotic systems. This 
also holds when one adds noise on the dynamics. A particularly useful set of tools is 
provided by ergodic theory and the thermodynamic formalism. In this approach one is 
interested in the statistical behavior of orbits, characterized by a suitable set of proba- 
bility measure. A natural choice are Gibbs measures in the sense of Sinai- Ruelle-Bowen 
[45]. In a forthcoming paper we indeed show that Gibbs measures arise naturally in 
BMS model. They come either from statistical inference principles where one tries to 
max;imize the statistical entropy given a set of fixed given quantities such as correla- 
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tions functions or mean firing rate (a prominent example of appUcation of this principle 
is given in [49]). They also arise when one wants to study the efi^ect of synaptic plas- 
ticity (learning, STDP) on the selection of orbits. In the context of BMS model one 
can show that Hebbian learning and STDP are related to a variational principle on the 
topological pressure, which is the analogon of free energy in statistical mechanics. 

5.6 The hmit dt -» 0. 

In the definition of the BMS model, one uses a somewhat rough approximation con- 
sisting in approximating the differential equation of the Integrate and Fire model with 
a Euler scheme, and discretizing time. A central question is: what did we lose by doing 
this, and is the model still relevant as a neural network model ? As mentioned in the 
introduction, this requires developments done elsewhere [14]. But we would like here 
to point out here a few remarks on this aspect. 

— Prom the "biological" point of view the Integrate and fire model with continuous 
time is already a rough approximation where the characteristic time for the neuron 
response is set to zero. One can actually distinguish (at least) 3 characteristic 
time scales in neuron dynamics descriptions based on difi^erential equations. The 
"microscopic time" dt corresponds somehow to the shortest time scale involved 
in the spike generation (e.g. microscopic mechanism leading to opening of ionic 
channels) . The "reaction time" Tr of the neuron corresponds to the time of raise and 
fall for the spike. If one focuses on spikes (and does not consider time averaging over 
sliding windows leading to the firing rate description) the last relevant time scale 
is the characteristic time T required for the neural network to reach a stationary 
regime. One expects to have dt « Tr « T . lu the IF model, however, the time 



reaction Tr is considered to be instantaneous (thus Tr < dt). This leads to delicate 
problems for the definition of the time of firing and requires the introduction of 
the "t~ notation". Using a discrete time approximation allows to circumvent this 
problem and corresponds somehow to pose dt = Tr = 1- 

One may reject this procedure a priori. Our philosophy is instead to extract as much 
results as possible from the discrete time spiking model and decide a posteriori what 
has been lost (or won). 

Prom the dynamical system point of view, the limit dt — » raises two problems. 
On one hand, the trajectories become continuous. Then one may have situations 
where the trajectory accumulates on S and where a small variation of the Wij's is 
not able to remove the intersection (as it is the case in th. 3). This type of situation 
is known in the field of genetic networks (see [23] and references therein). However, 
as mentioned in the paper, the situation is slightly different here, because of the 
neurons reset, leading to an infinite contraction of a domain onto a point. This 
effect really simplifies the dynamics study, and is still present in the continuous 
time case. However, this aspect would require careful investigations, not in the 
scope of the present work. 

The second problem is the use of a Euler scheme in the discretization. Using more 
elaborated schemes would complicate the analysis since the model would loose its 
convenient piecewise affine structure. We don't know what this would add. 
Finally, from a numerical point of view, softwares use discrete time. One aspect that 
interests us particularly is to know what are actually the computing capacities of 
the discrete time model compared to classical IF models and how much has been 
lost. 
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Fig. 1 Fig. la. Average value of the distance d{A, S) versus 7, C, for = 50. Fig. lb. A' = 100 
(in logj^Q scale). 



